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Saturation and Thermalization of the Magnetorotational Instability: 
Recurrent Channel Flows and Reconnections 

Takayoshi Sano^ and Shu-ichiro Inutsuka^'^ 
ABSTRACT 

The nonlinear evolution and the saturation mechanism of the magnetorotational 
instability (MRI) are investigated using three-dimensional resistive MHD simulations. 
A local shearing box is used for our numerical analysis and the simulations begin with 
a purely vertical magnetic field. We find that the magnetic stress in the nonlinear 
stage of the MRI is strongly fluctuating. The time evolution shows the quasi-periodic 
recurrence of spike-shape variations typically for a few orbits which correspond to the 
rapid amplification of the magnetic field by the nonlinear growth of a two-channel 
solution followed by the decay through magnetic reconnections. The increase rate of the 
total energy in the shearing box system is analytically related to the volume-averaged 
torque in the system. We find that at the saturated state this energy gain of the 
system is balanced with the increase of the thermal energy mostly due to the joule 
heating. The spike-shape time evolution is a general feature of the nonlinear evolution 
of the MRI in the disks threaded by vertical fields and can be seen if the effective 
magnetic Reynolds number is larger than about unity. 

Subject headings: accretion, accretion disks — diffusion — instabilities — MHD — 
turbulence 



1. INTRODUCTION 

The understanding of the origin of angular momentum transport has been required for 
the development of accretion disk theory. Balbus & Hawley (1991) have shown that weakly 
magnetized accretion disks are subject to a powerful local MHD instability. The discovery of 
this magnetorotational instability (MRI) brought about huge progress in the theoretical picture 
of accretion disks. Numerical simulations have revealed that the nonlinear evolution of the MRI 
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leads to MHD turbulence in which angular momentum is transported by the Maxwell stress (e.g., 
Hawley, Gammie, & Balbus 1995). Because the saturation amplitude of the stress determines 
the structure and evolution of the disks, the mechanism of the nonlinear saturation has a great 
importance for the theory. However the mechanism which regulates the saturation amplitude of 
the MHD turbulence is still unclear. 

The purpose of this letter is to examine what is happening at the saturated turbulent state 
in order to find a clue to the understanding of the saturation mechanism. A local shearing 
box approximation is used for our numerical calculations so the evolution of the instability can 
be followed for many orbital periods. According to previous numerical studies, the saturation 
amplitudes of the Maxwell stress take similar values for both local and global disk calculations 
(e.g., Stone et al. 1996). Thus we should find the essential processes of the saturation mechanism 
in the local shearing box simulations. 

The calculations include the effect of ohmic dissipation. Local 3D simulations show the 
nonlinear saturation of the MRI for both with and without the magnetic resistivity (Matsumoto 
&; Tajima 1995; Fleming, Stone, & Hawley 2000). In axisymmetric 2D ideal MHD simulations 
with initially uniform vertical field, the evolution of the MRI results in an exponentially growing 
solution (two-channel flow) with no nonlinear saturation (Hawley & Balbus 1992). When the effect 
of ohmic dissipation is efficient enough, however, the nonlinear saturation occurs even in 2D (Sano, 
Inutsuka, & Miyama 1998). These results suggest that the resistivity is an important quantity in 
determining the saturated state. In this letter we focus on the effect of ohmic dissipation in 3D 
calculations. This approach has an advantage of reducing the influence of numerical dissipation 
on the saturated state that is characterized by the physical (i.e., ohmic) dissipation process. 

2. NUMERICAL CALCULATIONS 

The resistive MHD equations are solved by using a finite-difference code developed by Sano, 
Inutsuka, Sz Miyama (1999). We use the local shearing box approximation described in detail 
by Hawley et al. (1995). The Cartesian coordinate {x, y, z) is defined around the fiducial point 
comoving with a fiducial angular velocity fi. Initially all the physical quantities are spatially 
uniform {p = po and P = Pq, where po and Pq are constant) except for the azimuthal velocity 
Vy{x) = ~qQx, where q = 3/2 for a Keplerian disk. In this letter, we focus on an initial magnetic 
field configuration with a uniform vertical field, = B^o, although the nonlinear evolution of the 
MRI are affected by the field geometries (e.g., Fleming et al. 2000). 

We choose normalizations with po = 1, Q, = 10~^ and the computational box has radial 
size Lx = 1, azimuthal size Ly = 4, and vertical size = 1. Most of the runs use a standard 
grid resolution of 32 x 128 x 32. The initial characteristic wavelength of the MRI is given by 
Amri = 27ri;A2o/^ for th^ ideal MHD case (Balbus &: Hawley 1991), where vazO = -Bzo/\/4vrpo is 
the Alfven speed. Its ratio to the vertical box size is assumed to be Xmri/ Lz = 0.32 for all models. 
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Then the initial field strength of the models is vazO = 5.0 x 10~^. Note that the net flux associated 
with the initial vertical field remains unchanged in the time evolution even if the resistivity is 
included. We assume the initial ratio of the gas and magnetic pressure is /3o = Po/iBzo/STr) = 10^, 
then the initial gas pressure is Pq = 1.25 x 10~^ and the sound speed is Cso = 4.56 x 10"'^ with 
7 = 5/3. The magnetic diffusivity r] is assumed to be constant and characterized by the magnetic 
Reynolds number Rcm- We define the magnetic Reynolds number as ReM = VL/rj = v\^Q/rjVL 
using typical velocity V = ua^o and typical length scale L = v^zo/^- 

3. LONG TIME EVOLUTIONS OF THE MRI 

The efficiency of angular momentum transport is given by the turbulent shear stress, 
Wxy = —BxBy/4:n + pVxSvy, where the first and second terms are the Maxwell and Reynolds stress, 
respectively, and Svy_= Vy + gfix. Figure 1 shows the time evolutions of volume-averaged Maxwell 
stress, {—BxBy/4:TT)B The top panel is the result of the case with ReM = 1 and the bottom is of 
the ideal MHD {rj = 0) case. 

The Maxwell stress dominates the Reynolds stress throughout the turbulent stage by a factor 
of about 6. The stress of the nonlinear stage is highly fluctuating for both cases shown in Figure 
1. The main feature of the time evolutions is the many spike-shape excursions in the stress. Power 
spectral density declines with frequency w approximately as at frequencies above a few orbital 
period, and almost flat at lower frequencies. Note that the break point corresponds to the typical 
timescale of the spike. 

The time averages of the Maxwell stress are shown in Figure 1 by dotted lines and they 
are 0.024 Pq and 0.013 Pq for the Rcm = 1 and r] = runs. Because the contribution of the 
spike-shape variations to the average is large, the understanding of the time evolution of these 
spikes is essential to estimate the saturation amplitude of the stress. 

The diffusion length l^m = r]/vAzO = 0.05 is larger than the grid scale Zgrid = 1/32 = 0.03 for 
the Rcm = 1 run, so that we can resolve the physical effect of the resistivity. In the ideal MHD 
run, on the other hand, the dissipation of the magnetic field is due to artificial numerical effects 
(Hawley et al. 1995). Although the resemblance of the time evolutions between the Rcm = 1 and 
r] = runs might suggest that the numerical dissipation works in a similar way to the physical 
resistivity, we use the results of the Rcm = 1 run in the following discussions. 

At the peak of each spike, a well-organized two-channel solution dominates in the entire 
computational box. Figure 2 shows images of the azimuthal component of the magnetic field, 
that is the dominant component at the saturated state. Top panel is the image at a time of 
the peak of a spike. The distribution is almost independent of x and y. The upper half of 



^The single brackets (/) imply a volume average of quantity /. We also use double brackets {(/)) to denote a time 
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fluid has a strong horizontal field with positive and negative By, and the lower half has the 
oppositely directed field with negative B^ and positive By. The perturbed velocity is nearly 
horizontal in the same direction as the fields but has the maximum at the neutral sheet of the 
fields. Because the gas pressure is still larger than the magnetic pressure by more than two orders 
of magnitude even at the nonlinear stage, the evolution is almost incompressible. Even at the 
maximum growth of the channel flow, the fluctuations of the density and pressure are very small, 
{5p^y/^/{p) ~ 7.2 X 10~3 and {5P^y/^/{P) ~ 1.5 x 10^^^ compared with the fluctuation of the 
magnetic pressure, ('5-Pmag)^''V(-Pmag> ~ 0.75. The exponentially growing channel solution is an 
exact solution for nonlinear incompressible MHD equations, thus the horizontal fields are amplified 
up to \Bx\ max ~ I8-B20 cind I -By I max ~ 30-6^0 by this nonlinear growth. The perturbed velocity is 
of the order of the Alfven speed at this time and is still subsonic; {\5v\) ~ 0.56(^;a) ~ 0.073(cs). 

The image at a time just after the peak is shown in the middle panel of Figure 2. The channel 
flow is known to be unstable for the parasitic instability (Goodman &: Xu 1994). Because the 
growth rate is proportional to the amplitude of the channel solution, this instability appears after 
the sufficient growth of the two-channel flow. The unstable modes have finite kx and ky and their 
wavelength 27r/(A:^ + ky)^^'^ must be longer than the vertical length of the channel L^. This kind 
of pattern occurs in the image of the spatial distribution of the magnetic field. The parasitic 
instability generates vertical motions and gives a chance for oppositely directed magnetic fields to 
approach each other. Then the dissipation of the magnetic field takes place efficiently through the 
magnetic reconnections. Finally the channel fiow disappears and a disorganized MHD turbulence 
lasts until the next channel fiow starts growing (bottom panel of Fig. 2). The magnetic energy 
is an order of magnitude smaller than the peak value and is almost equi-partitioned with the 
perturbed kinetic energy; {\Sv\) ~ 0.83(fA) ~ 0.053(cs). 



ENERGY BUDGET IN THE MHD TURBULENCE 



According to Hawley et al. (1995), we define the total energy within the shearing box as 

.,2 \ ^2- 



r = 



(1) 



where e means the specific internal energy and = —qQ'^x'^ is the tidal expansion of the effective 
potential. Using the evolution equation for the resistive MHD system, the time-derivative of the 
above equation gives 

^ = qnix dA ^pvjvy - ^) , (2) 

where dA is the surface element and the integral is taken over either of the radial boundary. Thus 
the increase rate of the total energy is proportional to the stress Wxy at the radial boundary. 
Balbus & Papaloizou (1999) have shown that the similar relation in the cylindrical coordinates 
holds for the global disk problem. Although the final expression of equation (§) does not explicitly 
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depend on the amount of resistivity, the joule heating is crucial for the transformation of the 
magnetic energy into the thermal energy (see below). 

Using the volume- averaged stress {wxy) instead of the surface- averaged value of the stress at 
the radial boundary {wxy)x, the volume average of the input energy (-Em) is approximately given 
by 

d(E- ) 

(Ein) = q'^ = {Wxy)x ^ 9^ {'^■xy) ■ (3) 

Because our numerical scheme solves the energy equation in terms of the total energy, our 
calculations perfectly satisfy equation (|2|). Furthermore, the density does not vary much so that 
the change in potential energy pcf) is negligible compared with the other terms in equation (jT|). 
Thus the input energy E^jn is almost identical to the increase of the sum of Sth = = P/{'^ ~ 1)) 
£^mag = ^VStt, and 

in — P't' /2- This is clearly shown in Figure 3 which depicts the time 
evolutions of {E\a) and (Etot) = {Eth) + {Emag) + {Ekm)- {Ein) and (-Etot) take such close values 
that their curves overlap each other. 

The time evolution of (Eth) is shown in Figure 3 by a dotted curve. Because (-Ein) is 
proportional to the Maxwell stress, its time evolution has spike-shape variations. As seen 
from the figure, the evolution of {Eth) also has many spikes with the similar amplitude to 
(-Ein). The time average of (-Eth) takes almost the same value as that of (-E^). The difference 
I ((-Eth)) — ((-Ein)) I /((-Ein)) is less than 10~^, when we calculate the time averages from 10 to 25 
orbits. Thus the energy gain of the system is finally transformed into the thermal energy. The 
thermal energy must be increasing throughout the nonlinear evolution unless some cooling process 
is included. The magnetic and kinetic energies, on the other hand, are saturated, thus the time 
averages of their rates of change are {{Eraa.g)) ~ ((-Ekin)) ~ 0- This means that the time and spatial 
average of the energy gain dT/dt should amount to the time and spatial average of the heating 
rate due to turbulent dissipation, provided the magnetic and kinetic energy are saturated. That is 

{{Eth)) « qH^xy)) ■ (4) 

Here we cannot remove the symbol (()) from this equation. In other words, this fluctuation- 
dissipation relation holds only in the coarse-grained description of the thermalization rate and the 
correlated fluctuation of magnetic field and velocity perturbations. Note that this equation does 
not hold in the 2D simulation with Rcm ^ 1 where the system does not show the saturation. 

Figure 4 shows the time evolutions of (-Ein), (-Eth)) and (-Emag) at a typical spike-shape 
variation. The joule heating rate (-Ejouie) = ??|V x B\'^/4it is also shown in this figure. The joule 
heating is about 80 percent of the thermal energy increase. (This fraction is almost 100 % for the 
Rgm = 0.3 run.) Throughout the nonlinear evolution, the magnetic energy is transformed into 
the thermal energy at a rate which over time matches (-Ein). Thus the saturation of the magnetic 
energy is achieved by a time-averaged balance between the energy gain by the growth of the MRI 
and the loss by the joule heating in the reconnection process. 

We find that the peak in (-Eth) is always delayed from (-Ein), and that the sign of (-Emag) 
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changes around the peak in (Eth)- The kinetic energy has a similar time evolution to -Emag; 
but the amplitude of the fluctuations in (E'kin) is negligible compared with that of the magnetic 
energy. The spike-shape variation of (i?in) starts with the onset of a channel flow. The magnetic 
energy increases due to the growth of the MRI and {E^g,g) takes a positive peak value at the 
maximum growth of the two-channel solution, which corresponds to the peak in {Em) (top panel 
of Fig. 2). When the channel flow is destroyed by the parasitic instability, the large scale magnetic 
reconnections become efficient. Then the sign of (Emag) changes from positive to negative and the 
thermal energy increases rapidly. When the (-Emag) takes a negative peak, (-Eth) is nearly at the 
top of the peak (middle panel of Fig. 2) . 



5. DEPENDENCE OF THE SATURATION CHARACTER ON THE MAGNETIC 

REYNOLDS NUMBER 

We find that two-channel solutions appear repeatedly in nonlinear turbulent stage. The 
behavior similar to this appearance of a channel fiow can be seen in 2D simulations, which show 
a two-channel fiow appears at the end of the calculations even though the initial value of Amri 

is smaller than the box size. This inverse cascade is characteristic of both 2D and 3D nonlinear 
evolutions of the MRI in the cases Rcm ^ 1- The difference between 2D and 3D simulations is 
that in 3D the channel flow is unstable to non-axisymmetric modes of the parasitic instability. 
Then the question is why the two-channel solution appear in nonlinear stage of the MRI. 

The key quantity is the characteristic length of the instability expected from the linear 
analysis. The most unstable wavelength of the MRI is given by Amri ~ va/^ oc va for less 
resistive cases Rcm ^ 1- If the field is amplified by the growth of the MRI, then Amri shifts to 
the longer wavelength. Therefore it can reach the vertical box size, and the longest wavelength 
mode, a two-channel solution, appears finally in nonlinear evolution. This is the reason for the 
inverse cascade. For Rcm ^ 1, on the other hand, the characteristic wavelength is given by 
Amri ~ v/'^^a (Sano &: Miyama 1999). Thus the wavelength Amri becomes shorter when 

the field is amplified. In fact, 2D simulations with Rcm ^ 1 show that the saturated state is a 
disorganized MHD turbulence and no emergence of a two-channel fiow (Sano et al. 1998). We 
find this is true even for 3D cases. For the case with Rcm = 0.3, the nonlinear evolution shows 
no growth of a two-channel flow. The amplitude of the time variability of the Maxwell stress is 
very small in this case. Due to the lack of the nonlinear growth of the channel flow, the horizontal 
component of the magnetic energy is not amplified very much and remains only 2 times larger 
than the vertical, while this ratio {{By)) / {{BD) is an order of magnitude larger for Rcm ^ 1 cases. 
The energy balance at the saturated state is roughly given by {{Em)) ^ ((-E'th)) ^ ((-Ejouie)) same as 
in the cases with Rsm ^ 1- Because the diffusion time is comparable to the growth time of the 
MRI in this case, the magnetic diffusion would be an important process for the joule heating as 
well as the recurrent and rapid magnetic reconnections. The saturation amplitude is two orders of 
magnitude smaller than the models with Rcm ^ 1) so that the angular momentum transport by 



- 7- 



the MRI is not efficient wlien Rbm ^ 1- 

We can categorize tlie saturated states of tlie MRI into two types wfien tfie disks are tfireaded 
by uniform vertical fields initially. The first type involves the recurrence of the channel flow 
and large scale reconnections (the Rbm = 1 and ideal MHD runs). In the disorganized MHD 
turbulence of the second type, the magnetic diffusion and reconnections prevail without the 
significant growth of the channel flow (the Rcm = 0.3 run). The angular momentum transport by 
the magnetic stress is efficient only for the first type. The difference of these two types is caused 
by the different dependence of the characteristic length of the MRI, Amri, on the magnetic field 
strength. Thus the magnetic Reynolds number defined as Rsm = Va/v^ can characterize the 
linear and nonlinear features of the MRI and the critical value is Rcm ~ 1- 

Wc thank James Stone and Neal Turner for useful discussions. Computations were carried 
out on the VPP5000 at the National Astronomical Observatory of Japan and the VPP700 at the 
Subaru Telescope, NAOJ. 
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Fig. 1. — Time evolution of volume- averaged Maxwell stress for the model with Rcm = 1 (top) 
and 7] = {bottom). Time is given in the rotation time trot = Dotted lines show their 

time-averaged values done from 10 to 100 orbits. The grid resolution of these runs is 32 x 128 x 32. 
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Fig. 2. — Slices in the x-z plane at y = —2 and in the y-z plane at x = —0.5 of the azimuthal 
component of the magnetic fields By/BzO {colors) and magnetic field vectors (arrows) in the 
Rcm = 1 run with grid resolution 64 x 256 x 64. From top to bottom the time of the image is 
t/trot = 21.5, 21.7, and 22.4. 
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Fig. 3. — Time evolution of volume-averaged time derivative of input energy (.Ein), total energy 
(E'tot) = {Eth) + (-Emag) + {Ekin), and thermal energy (E^th) of the Rbm = 1 run with grid 
resolution 64 x 256 x 64. The time derivatives of the energies are given in units of £^tho/irot 
where E'tho = -fb/(7 — 1) is the initial thermal energy. 
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Fig. 4. — Time evolution of volume-averaged time derivative of input energy (£'in), thermal energy 
(Eth), and magnetic energy (-Emag) and volume-averaged joule heating rate (£^jouie) at a spike- 
shape variation in Reu = 1 run with grid resolution 64 x 256 x 64. Arrows denote the times of the 
snapshots shown in Figure 2. 



